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We report the first simulations of the Faraday instability using the full three-dimensional 
Navier-Stokes equations in domains much larger than the characteristic wavelength of the 
pattern. We use a massively parallel code based on a hybrid Front-Tracking/Level-set 
algorithm for Lagrangian tracking of arbitrarily deformable phase interfaces. Simulations 
performed in square and cylindrical domains yield complex patterns. In particular, a 
superlattice-like pattern similar to those of [Douady & Fauve, Europhys. Lett. 6, 221-226 
(1988); Douady, J. Fluid Mech. 221 , 383-409 (1990)] is observed. The pattern consists 
of the superposition of two square super lattices. We conjecture that such patterns are 
widespread if the square container is large compared to the critical wavelength. In 
the cylinder, pentagonal cells near the outer wall allow a square-wave pattern to be 
accommodated in the center. 


1. Introduction 


In the Faraday (1831) problem, the interface between two superposed fluid layers, sub¬ 


jected to periodic vertical vibration of sufficient amplitude, forms sustained standing wave 
patterns, called Faraday waves. Many universal dynamical-systems phenomena were first 
discovered through the Faraday experiment. In particular, a number of experimental in¬ 
vestigations have focused on the variety of patterns formed in square (e.g. |Douady ~fc 


& Gollub, 1985; Das & Hopfingerl 2008; Batson, Zoueshtiagh & Narayanan, 2013; Ra- 

jchenbach, Clamond & Leroux 

2013) containers, which have led to important theoretical 

developments (e.g. Meron 1987 

Silber & Knobloch 1989 Crawford, Knobloch & Riecke, 


tern formation. 


Recently, Perinet, Juric & Tuckerman (2009 2012) performed the first three-dimensional 


direct numerical simulation of Faraday waves in horizontally periodic domains contain¬ 
ing a few basic cells. Here, we describe numerical simulations of Faraday waves in much 


f Email address for correspondence: laurette@pmmh.espci.fr 




























































2 L. Kahouadji, N. Perinet, L.S. Tuckerman, S. Shin, J. Chergui and D. Juric 


larger domains which have yielded more complicated patterns. For this purpose, we have 
used the free-surface code developed by |Shin, Chergui fe Juric (2014) for simulations 
of two-phase flows on massively parallel computer architectures. In a square domain 
(i.e. with square horizontal cross-section), we have obtained patterns with two different 
length scales, very similar to those reported in experiments by Douady & Fauve (1988) 


and Douady (1990). Patterns of this type were later termed superlattices by Kudrolli, 


Pier & Gollub (1998) and Arbell & Fineberg ( 2002 ), who observed a large variety of such 
patterns in their experiments, but by using a temporal forcing with two frequencies. The 
spectral analysis of our patterns, produced by single-frequency forcing, shows that they 
result primarily from the superposition of two superlattices constructed from very similar 
wavelengths. We have also simulated the same conditions in a cylindrical container. The 
resulting interface contains squares and pentagons. 

The paper is organized as follows. First, we briefly present the problem and explain 
the methods involved in the code. We then present the numerical linear stability analysis 
that is validated by comparison with the theoretical treatment of |Kumar fe; Tuckerman| 
(1994). We demonstrate the robustness of the obtained superlattice-like pattern by 
testing several boundary conditions. This pattern is described and analysed spectrally. 
Then, we show the results of simulations with the same conditions but in a cylindrical 
domain. 


2. Problem formulation, governing equations and numerical scheme 

The governing equations for an incompressible two-phase flow can be expressed by a 
single field formulation: 

p^ + u-Vuj = -VP + pG + V •(Vu + Vu T ) +F, V • u = 0 (2.1) 

where u is the velocity, P is the pressure, p is the density, p is the dynamic viscosity, G 
represents the homogeneous volume forces, and F is the local surface tension force at the 
interface. Here, G contains the gravitational term supplemented by a time-dependent 
force accounting for the vibrations of the domain: 

G = — (g + 7 cos(u;t)) e z ( 2 . 2 ) 


where g is the gravitational acceleration, e z is the upward vertical unit vector, 7 is the 
amplitude of the inertial forcing and uo = 27r/T its frequency. 

The code essentially consists of two modules: one, which solves the three-dimensional 
incompressible Navier-Stokes equations, and the other, which treats the free surface using 
a parallel Lagrangian Tracking method and is only active if the flow is a two-phase flow. 
Material properties such as density or viscosity are defined in the entire domain: 


p(x, t) * Pi 4 - (p 2 - pi ) /(x, t) 

m(x, 0 = //, • (/i 2 - Mi) ^( x > t)- 


(2-3) 


The indicator function, /, is a numerical Heaviside function, ideally zero in one phase 
and one in the other phase. I is resolved with a sharp but smooth transition across 3 to 
4 grid cells and is generated using a vector distance function computed directly from the 


tracked interface (Shin & Juric (2009a)). 

The fluid variables u and P are calculated by a projection method ( Chorin| |1968 ). 
The temporal scheme is first order, with implicit time integration used for the viscous 
terms. For spatial discretization we use the staggered-mesh marker-in-cell (MAC) method 
(Harlow & Welch 1965) on a uniform finite-difference grid with second-order essentially 
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non-oscillatory (ENO) advection (Shu & Osher 1989). The pressure and distance func¬ 
tion are located at cell centers while the x, y and z components of velocity are located 
at the faces. All spatial derivatives are approximated by standard second-order centered 
differences. The treatment of the free surface uses a hybrid Front-Tracking/Level-Set 
technique which defines the interface both by a discontinuous density field on the Eule- 
rian grid as well as by triangles on the Lagrangian interface mesh. 

The surface tension F is implemented by the hybrid formulation 


F = Kh = 


F l -N 

N-N 


(2.4) 


where a is the surface tension coefficient assumed to be constant and kh is twice the 
mean interface curvature field calculated on the Eulerian grid, with 


F L= K f n f s f ( x - x /) ds, N = / n f 6f (x - X/ ) ds 

Jr(t) Jr(t) 


(2.5) 


Here, xj is a parameterization of the time-dependent interface, T(t), and Sf (x — xj) is 
a Dirac distribution that is non-zero only where x = xj; Uf stands for the unit normal 
vector to the interface and ds is the length of an interface element; Kf is twice the mean 
interface curvature obtained on the Lagrangian interface. The geometric information, 
unit normal, n/, and interface element length, ds in N are computed directly from the 
Lagrangian interface and then distributed onto the Eulerian grid using the discrete delta 
function and the immersed boundary method of Peskin (1977). A detailed description 


of the procedure for calculating F, N and I can be found in^Shin (2007). 


The Lagrangian interface is advected by integrating dx //dt = V with a second-order 
Runge-Kutta method where the interface velocity, V, is interpolated from the Eulerian 
velocity. 

The time step At is chosen at each iteration in order to satisfy a criterion based on 


{AtcpL, Ati n t, At v i s , 

^Gap} 

which ensures stability of the calculations. These bounds are defined by: 


A Xj 


A^cfl = min ( min 

j ydomain y ujj 

At vis = min f —, — j 

VM2 Mi/ 6 


Atint = min ( min 


( Ax j \ 


At, 


cap 


j v r (*) v IV|| / 

1 f {pl + P2A x\ 

2 V 7TCT 


Ln 


1/2 


( 2 . 6 ) 


(2.7) 


where Ax™in = 


(A Xj). With the periodic volume force term G of (2.2), a supple- 
27r/(50cd) is required. 


mentary upper bound A t u 

Fortran 2003 allows the definition of a set of dynamically allocated derived types 
and generic procedures associated with the grids, scalar and vector fields, operators as 
well as the various solvers used in the Navier-Stokes and Lagrangian tracking modules. 
The parallelization of the code is based on algebraic domain decomposition, where the 
velocity field is solved by a parallel generalized minimum residual (GMRES) method for 
the implicit viscous terms and the pressure by a parallel multigrid method motivated by 
the algorithm of Kwak & Lee (2004). Communication across process threads is handled 
by message passing interface (MPI) procedures. 

Finally, the code also contains a module for the definition of immersed solid objects 
and their interaction with the flow, which we have used to simulate Faraday waves in a 
cylindrical container. A Navier-slip dynamic contact line model is implemented where 
hysteresis is accounted for by fixing the contact angle limits to prescribed advancing or 
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Figure 1. a) Evolution of the total kinetic energy for uo/2n = 30 Hz and various accelerations 
in a small domain. The fact that the peaks for a given acceleration lie along a line shows that 
the energy growth is very close to exponential, b) Evolution of the various bounds on the time 
step during the course of the simulation. The smallest, At C a P , is used in the simulations. 


uj/2tv 

Hz 

Ac 

mm 

Floquet 

7c/s 

DNS 

7c/ g 

relative error 
Floquet/DNS 

30 

11.3 

0.733 

0.713 

2.73 % 

60 

5.9 

2.65 

2.61 

1.50 % 

90 

4.3 

5.32 

5.3 

0.38 % 


Table 1 . Comparison of Floquet and DNS instability thresholds for various frequencies. 


receding angles as in Shin & Juric (2009b). Further details are available in Shin, Chergui 


& Juric (2014). 


3. Parameters and linear stability 

We consider a layer of silicone oil of thickness h = 14.5 mm with density pi = 965 
kg m -3 and dynamic viscosity pi = 0.02 kg m _ 1 s _1 . This layer is overlaid with 29.5 
mm of air of density p 2 = 1.205 kg m -3 and dynamic viscosity p 2 = 1.82 x 10 -5 kg 
m - 1 s -1 . The surface tension at the interface between the oil and the air is a — 0.02 
kg s -2 . This choice of parameters was originally motivated by experiments on Faraday 
waves in a rotating cylindrical container (Clement, Pucci & Couder, unpublished), which 
will be the subject of a future investigation. 

Before simulating the patterns, we computed the critical acceleration y c from direct 
numerical simulations (DNS) and compared it with that of the Floquet linear stability 
analysis of Kumar & Tuckerman (1994). These simulations are conducted in a small 
periodic domain of length A c and width A c / 2 , where the critical wavelength is A c = 
27r/k c = 11.3 mm with k c the critical wavenumber. We compute the growth rate of the 
total kinetic energy of the system (figure [lj a)) after the growth becomes exponential. 
Near the threshold y c , the growth rate varies linearly with 7 , so that y c can be deduced 
by linear interpolation. Table [l] compares the Floquet predictions of y c to the results 
from DNS for three vibration frequencies. The discrepancy is under 3% in all three cases. 
The ratio of gravitational to capillary forces is measured by the Bond number, whose 
value here is gAp/(ak 2 ) = 1.5, indicating that both forces are at work. 
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Figure 2. Calculation domain of size L x L x L/3 (L = 132 mm), divided into 8x8x4 
subdomains, a) Interface profile, b) an immersed solid cylinder of diameter D — L and height 
h — L/3 surrounding the flow. The flow outside of the cylinder is not computed. For both 
(a) and (b) the Cartesian resolution per subdomain is 64 3 , which gives a global resolution of 
512 x 512 x 256. 



Figure 3. Snapshots of the top view of the interface colored by the vertical velocity for (a) 
no-slip, (b) free-slip, and (c) periodic boundary conditions. 


Our large-scale simulations are carried out with a vibration frequency of u;/27r = 30 
Hz, i.e. a period of T « 0.033 s, and an amplitude of 7 = g = 1 . 367 c . Figure 0 b) 
tracks the timestep bounds (2.7) during the simulation. It can be seen that At cap is 
the limiting timestep throughout the simulation and is on the order of 5 x 10 -5 s. The 
domain is of size L x L x L/ 3, where L = 132 mm = 11.7 A c . This domain is subdivided 
into 8 x 8 x 4 = 256 subdomains whose traces can be seen in figure [2|a), each of which 
contains 64 3 gridpoints, leading to an overall resolution of 512 x 512 x 256 gridpoints 
on a Cartesian mesh. The critical wavelength is thus resolved by 44 = 512A C /L grid 
cells. Each subdomain is assigned to a process thread. Our initial condition has zero 
velocity and random noise on the interface of order 10 -2 mm. The pattern emerges 
after approximately 40 forcing periods, increases exponentially, and saturates at about 
90 forcing periods, requiring about 15 days of computation time on 256 threads of an 
IBM x3750-M4. The cylindrical container has been modeled by an impermeable solid 
object as shown in figure [ 2 jb). 


4. Results and discussion 

4.1. Visualisations 

The top and bottom of the domain are taken to be rigid; we impose no-slip boundary 
conditions on the velocity field. Three types of boundary conditions have been tested 
on the lateral walls, (a) no-slip, (b) free-slip, and (c) periodic. Although our parameters 
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Figure 4. Interface profiles at times separated by intervals of T/2, where T is the forcing period 
and 2 T is the subharmonic response period. In (a) and (c), the small squares consist of wells 
surrounded by ridges with peaks at each corner, while in (b), low central peaks are surrounded 
by four higher peaks. In (a), a weak diagonal “bridge” connects the large squares on the bottom 
left and top right; in (c), the bridge links the top left and bottom right. 



Figure 5. Velocity profiles on the plane y — 3L/4, corresponding to each of the three instants 
shown in figure^] The vectors are colored according to their vertical component. 


were originally chosen for a different purpose, we find that the patterns that we observe 
are strikingly similar to those observed in earlier experiments by |Douady fe Fauve (1988), 
figure 1(b) and Douady (1990), figures 10(b), 11(b) and 15(b), with different parameters: 
water, with pi = 1000 kg m -3 , pi = 9 x 10 -4 kg m -1 s -1 , and a = 0.03 kg s -2 , in a 
container of dimensions 80 mm x 80 mm x 5 mm oscillated with / = 37 Hz or / = 70 Hz. 

For all three boundary conditions, we observe superlattice-like patterns shown in figure 
[3j The interface in the x-y plane consists of four large square sub-blocks, each composed 
of smaller squares whose sides have approximate length A c . In each case, the overall 
pattern has symmetry , i.e., it is invariant under reflections through the two diagonals. 
Each of the four sub-blocks is in phase opposition with its two adjacent neighbors. The 
patterns in (a) and (c) have periodicity length L while (b) has a period of 2 L. Although 
one would expect the periodic case (c) to be the simplest to analyze, we observe that the 
blocks in the free-slip case (b) are more homogeneous than those in cases (a) and (c); we 
will therefore focus on this case. 

Figure [4] displays the temporal evolution of the interface at saturation. (A movie of 
this evolution is available as supplementary material to this article.) After one period 
of forcing oscillation T, troughs are replaced by crests (see figure [4ja) and (c)) as a 
consequence of the subharmonic behavior of the interface. This subharmonic behavior is 
also displayed by the large sub-blocks, since the appearance of each block is transformed 
into that of its two adjacent neighbors after time T and then returns to its initial shape 
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a) 




Figure 6. Spatial spectrum of a snapshot of the interface in the doubled domain. Circles C\ 
(dark blue) and C 3 (dark brown) have slightly different radii. C\ contains the fundamental mode 
£23,1 and its images under reflection and rotation (dark blue dots surrounded by circles). C 3 
contains the highest third-harmonic mode £23,3 and its images (dark brown dots surrounded by 
squares). Only these modes and, to a lesser extent, the weaker second harmonic £ 46,0 (orange) 
have significant amplitude, a) Restriction to [—25, 25] 2 highlights the square symmetry of the 
pattern, b) Restriction to the first quadrant [0,50] 2 includes £ 46 , 0 - 


after time 2 T has passed. The pattern is thus invariant under the combined operations 
of rotation by n/2 and time-translation by T, a spatio-temporal symmetry. 

Figure [5] shows the interface and the velocity on a vertical slice at the same three 
instants. The contrasting behavior of the velocity in the left and right halves of the 
slice is displayed. The interface of figure [5ja) has a maximum (minimum) at the left 
(right) boundary and vice versa for (c), illustrating again that the field is not periodic 
in a domain of length L. 


4.2. Fourier spectra 

For free-slip boundary conditions, the pattern is periodic when it is doubled horizontally 
by reflection at the boundaries. The doubled pattern is invariant under the symmetry 
group D 4 of the square, consisting of reflections through the x and y coordinate axes 
and rotations of multiples of 7r/2. We decompose the doubled interface height profile 
£(x, 2 /,£) into spatial Fourier modes, and then into spatio-temporal modes as follows: 


C(x,y,t) 


E Cj,mWe ikl -"*- x 

l,rn= — oo 


E 

l r m =—00 


E f inQt /2 


lk^ jT) 


(4.1) 


where k^ m = ( 77 , ^-p), x = and = |k/ ?m |. The square symmetry of the 

doubled interface implies that the eight modes in an octet C±i,±m and £±™,±z ad have 
the same amplitude. Figure ^ shows two such octets, with representative elements £ 23,1 
and £ 23 , 3 ? located on circles which we will call C\ and C3, respectively, i.e. with l = 23 
and m = 1 or m = 3. 

The circle C\ has radius fc c , the critical wavenumber. The dominant wavevector octet 
of our pattern, represented by k 2 3,i, is in accordance with the experimental observations 
of Douady & Fauve (1988), who scanned over the forcing frequency in order to vary the 
selected wavenumbers and patterns. They noted that the range of existence was greatest 
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for pairs (Z,ra) with a small ratio m/l and was maximum for m = 1 , as in our case, i.e. 
for wavevectors which are almost parallel to the walls. 

The modes on C 3 can be seen as the result of cubic interactions between modes of C \: 


*23 ,3 


= k 23 ,i + k 23 ,i + k_ 


-23,1 


(4.2) 


Modes on C 3 are fed continually and damped slowly, if at all, since ^ 23,3 — k c \/k c = 0.75%. 
Generalizing ( |4.2| ) to any m l shows that such cubic interactions lead to k/ ? 3 m , with 
|^z, 3 m ~ k c \/k c ~ 4 (m/l) 2 . Due to its proximity to the critical circle, this mode is weakly 
damped and should often be present when Faraday waves are excited in a square domain 
which is much larger than the critical wavelength. A second harmonic of C \, arising from 
the quadratic interaction 

k46,0 = k23,l + k 2 3,-l (4.3) 

is seen in figure [6jb) along with its image under 7 t/2 rotation. These modes occur in a 
quartet rather than an octet because of their reflection symmetry in the coordinate axes 
and are of smaller amplitude than the modes on C\ or C3. The experimental techniques 


of Douady & Fauve (1988) did not give them access to amplitudes of Fourier modes other 


than the dominant one. 


According to the terminology of |Simonelli fe Gollub (1989), the four co mponents 


C±i,±m comprise one pure mode, while C±m,±z comprise another pure mode. Douady 


& Fauve (1988) observed mixed modes (which they called symmetric) for small m/l, 


and pure modes (which they called dissymmetric) for larger m/l. Pure modes may be 
combined to form two types of mixed modes, the in-phase mode in which all of the com¬ 
ponents have the same sign (as in our figure [6| , and the out-of-phase mode in which 


the components of the two pure modes are of opposite sign. Silber & Knobloch (1989) 


showed that these two types of mixed modes are equivalent if l and m have opposite par¬ 
ities. Crawford (1991) showed that this was also true for l and m with the same parity 
(as in our case) by invoking the concept of hidden symmetries. In a square container 
with Neumann boundary conditions, the doubled pattern formed by reflecting at each 
boundary, as we have done above, satisfies periodic boundary conditions, whose inherited 
(hidden) translation invariances are incorporated into the normal form. 

For each of C\ or C3, the resulting pattern has a superlattice structure, in which squares 
of different sizes coexist, as shown in figure [TJa) and (b). Although shown in domains of 
length L, the periodicity length is 2 L. The smallest squares have length scale 27r/&23,m 
for m = 1,3. The patterns formed by modes on C 3 have an additional intermediate length 


scale 2L/3. Close examination of figure 
almost, but not exactly identical. Figure 


7T>) shows that adjacent medium squares are 
7[c) shows a superposition of the two patterns 


of figure^a) and (b), weighted by the amplitudes of their respective Fourier components. 
During most of the time, the interface resembles this superposition; compare with figure 
[4j When the amplitudes of modes on C\ and C 3 approach zero, which occurs twice during 
one subharmonic period (see figure J 8 ^a) and figure |4jb)), it is necessary to take modes 
like £ 46,0 into account. 

The temporal evolution of £ 23,15 £ 23,3 and £ 46,0 is displayed in figure Ji](a) and the 
spatio-temporal spectrum is shown in figure [ 8 jb). The mode amplitudes ( 23,1 and £ 23,3 
are related by a multiplicative constant and contain only subharmonic temporal compo¬ 
nents (i.e. n odd in ( |4.1| )). The temporal component n = 1 corresponding to D = uj/2 
dominates the temporal evolution. Both £ 23,1 and £ 23,3 evolve in phase opposition, cross¬ 
ing zero at the same time and yielding patterns such as that in figure [4^b). In contrast, 
£ 46,0 is harmonic and dominated by two temporal components n = 0 and n = 2, as ex- 




















Figure 7. Patterns constructed from Fourier modes (a) in C 1 , (b) in C 3, and (c) a superpo¬ 
sition of the two. Lighter (darker) zones correspond to the peaks (troughs) of the interface. 
Superlattices, consisting of smaller and larger squares, are observed in all of the plots, Smallest 
boxes are approximately of size A c . b) Medium squares of size 2L/3 are present (long-dashed 
box), c) Reconstructed pattern, made by combining patterns in (a) and (b) weighted by the 
corresponding Fourier coefficients £ 23,1 and £ 23 , 3 - Compare with figure Qc). 



Figure 8. a) Temporal evolution of the three main spatial modes. Modes £23,1 (blue curve 
and circles) and £ 23,3 (red curve and + symbols) are subharmonic: their period of oscillation is 
twice that of the forcing. Mode £ 46,0 (black curve and x symbols) is harmonic, b) Temporal 
spectrum of the main modes. (23,1, C23,3 contain only odd temporal harmonics while C46,o has 
mainly even harmonics. Each spatial mode has a finite n — 0 component and hence a non-zero 
temporal mean. 


pected from the quadratic resonance of £ 23,1 with itself or with ( 23 , 3 . Each spatial mode 
has a finite n = 0 mean component; the offset of £46,0 is especially noticeable. 

4.3. Simulation in a cylindrical container 

We performed simulations under the same conditions, but in a cylindrical container of 
diameter D = L and height h = L/ 3 (see figure [2|b)). The lateral boundary condi¬ 
tions are free-slip and the advancing and receding contact angles are fixed at 100° and 
60° respectively. Figure [ 9 ] shows snapshots of the interface every Tj 2 . (A movie of 
this evolution is available as supplementary material to this article.) The pattern has 
the spatio-temporal symmetry of invariance under combined time-evolution by T/2 and 
rotation by 7r/2, while each instantaneous pattern is reflection-symmetric about both di¬ 
agonal axes. At the center, we observe squares whose characteristic wavelength remains 
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Figure 9. Interface profile in the cylinder at time intervals of T/2. The central part of the 
domain is still occupied by squares, but five-sided cells are present near the boundaries, allowing 
the pattern to fit inside a circle. 


close to A c of Table [l] for / = 30 Hz. Closer to the boundary, five-sided cells can be seen 
as the pattern organizes itself to fit in the circular domain. Aside from this, the pattern 
is relatively insensitive to the circular shape of the domain because L/\ c = 11.7 1, 

in contrast with small-radius experiments and theory. In most of these (e.g. Ciliberto 

|Ba~n, 


fc Gollub[ |1985| Crawford, Knobloch fc Riecke) |1990| |Das Sz Hopfinger] |2008 
Zoueshtiagh & Narayanan, 2013), the interface takes the form of Bessel functions, includ¬ 


ing axisymmetric modes, while Rajchenbach, Clamond & Leroux (2013) report intriguing 
patterns of the form of stars and pentagons. 


5. Conclusion 

We have presented the first three-dimensional numerical results of the Faraday insta¬ 
bility in domains much larger than the critical wavelength in both square and cylindrical 
containers. In the square domain, the interface displays patterns that are very similar 
to those found in the work of Douady & Fauve (1988) and Douady (1990). By means 


of a spatial spectral decomposition, we have described these patterns quantitatively and 
have found that they have a complex superlattice-like structure due to resonances be¬ 
tween modes of very similar wavelengths. We conjecture that this pattern arises from 
two effects: (i) In a square container which is large compared to the critical wavelength, 
Faraday wave patterns tend to be mixed modes whose wavevectors are almost parallel to 
the boundaries of the box, as stated by Douady & Fauve (1988]). (ii) In this case, cubic 
nonlinearities generate another wavelength which is very close to the dominant one. The 
combination of the two effects leads to the superlattice pattern that we observe in our 
simulations and we infer that such patterns should be observed in a wide variety of sys¬ 
tems. The present work demonstrates that it is now possible to numerically compute all 
the fields of Faraday waves in domains sufficiently large to produce complex patterns such 
as super lattices, quasicrystalline patterns or oscillons. The hydrodynamic mechanisms 
responsible for the formation of these patterns can then be explored numerically. 


Acknowledgements 

We thank S. Douady, S. Fauve, and G. Pucci for helpful discussions. This work was 
performed using high performance computing resources provided by the Institut du De- 
veloppement et des Ressources en Informatique Scientifique (IDRIS) of the Centre Na¬ 
tional de la Recherche Scientifique (CNRS), coordinated by GENCI (Grand Equipement 
National de Calcul Intensif). N. Perinet was supported by a grant from CONICYT - 
FONDECYT/postdoctoral research project 3140522. L.S.T. acknowledges support from 







































11 


Numerical simulation of super-square patterns in Faraday waves 

the Agence Nationale de la Recherche (ANR) for the TRANSFLOW project. This re¬ 
search was supported by the Basic Science Research Program through the National Re¬ 
search Foundation of Korea (NRF) funded by the Ministry of Science, ICT and future 
planning (NRF-2014R1A2A1A11051346). 


References 

Arbell, H. & Fineberg, J. 2002 Pattern formation in two-frequency forced parametric 
waves. Phys. Rev. E 65, 036224. 

Batson, W., Zoueshtiagh, F. & Narayanan, R. 2013 The Faraday threshold in 
small cylinders and the sidewall non-ideality. J. Fluid Mech. 729, 496-523. 

Chorin, A. J. 1968 Numerical simulation of the Navier-Stokes equations. Math. Corn- 
put. 22, 745-762. 

Ciliberto, S. & Gollub, J. P. 1985 Chaotic mode competition in parametrically 
forced surface waves J. Fluid Mech. 158, 381-398. 

Clement, F., Pucci, G. & Couder, Y. Instability de Faraday en rotation, unpub¬ 
lished. 

Crawford, J. D., Knobloch, E., Riecke, H. 1990 Period-doubling mode interac¬ 
tions with circular symmetry Physica D 44, 340-396. 

Crawford, J. D. 1991 Normal forms for driven surface waves: Boundary conditions, 
symmetry, and genericity Physica D 52, 429-457. 

Das, S. P. & Hopfinger, E. J. 2008 Parametrically forced gravity waves in a circular 
cylinder and finite-time singularity J. Fluid Mech. 599, 205-228. 

Douady, S. & Fauve, S. 1988 Pattern selection in Faraday instability. Euro- 
phys. Lett. 6, 221-226. 

Douady, S. 1990 Experimental study of the Faraday instability J. Fluid Mech. 221, 
383-409. 

Faraday, M. 1831 On a peculiar class of acoustical figures; and on certain forms assumed 
by groups of particles upon vibrating elastic surfaces. Philos. Trans. R. Soc. London 
121, 299-340. 

Gomes, M.G.M. & Stewart, I. 1994 Steady PDEs on generalized rectangles: a change 
of genericity in mode interactions. Nonlinearity 7, 253-272. 

Harlow, F. H. & Welch, J. E. 1965 Numerical calculation of time dependent viscous 
incompressible flow of fluid with free surface. Phys. Fluids. 8, 2182. 

Kudrolli, A., Pier, B. & Gollub, J.P. 1998 Superlattice patterns in surface waves. 
Physica D 123, 99-111. See arXiv:chao-dyn/9803016 for more legible versions of fig¬ 
ures. 

Kumar, K. & Tuckerman, L. S. 1994 Parametric instability of the interface between 
two fluids. J. Fluid Mech. 279, 49-68. 

Kwak, D. Y. & Lee, J. S. 2004 Multigrid algorithm for the cell-centered finite- 
difference method II: Discontinuous coefficient case. Numer. Meth. Part. Dif¬ 
fer. Equ. 20, 723-741. 

Meron, E. 1987 Parametric excitation of multimode dissipative systems. 
Phys. Rev. A 35, 4892-4895. 

Perinet, N., Juric, D. & Tuckerman, L. S. 2009 Numerical simulation of Faraday 
waves. J. Fluid Mech. 635, 1-26. 

Perinet, N., Juric, D. & Tuckerman, L. S. 2012 Alternating hexagonal and striped 
patterns in Faraday waves. Phys. Rev. Lett. 109, 164501. 

Peskin, C. S. 1977 Numerical analysis of blood flow in the heart. J. Comput. Phys. 25, 
220-252. 


12 L. Kahouadji, N. Perinet, L.S. Tuckerman, S. Shin, J. Chergui and D. Juric 

Rajchenbach, J., Clamond, D., Leroux, A. 2013 Observation of star-shaped surface 
gravity waves. Phys. Rev. Lett. 110, 094502. 

Shin, S. 2007 Computation of the curvature field in numerical simulation of multiphase 
flow. J. Comput. Phys. 222, 872-878. 

Shin, S. & Juric, D. 2009 A hybrid interface method for three-dimensional multiphase 
flows based on front-tracking and level set techniques. Int. J. Num. Meth. Fluids 60, 
753-778. 

Shin, S. & Juric, D. 2009 Simulation of droplet impact on a solid surface using the 
level contour reconstruction method. J. Mech. Sci. Technol. 23, 2434-2443. 

Shin, S. Chergui, J. & Juric, D. 2014 A solver for massively parallel direct nu¬ 
merical simulation of three-dimensional multiphase flows. Comput. Fluids , submitted. 
arXiv: 1410.8568 physics.flu-dyn]. 

Shu, C. W. & Osher, S. 1989 Efficient implementation of essentially non-oscillatory 
shock capturing schemes, II. J. Comput. Phys. 83, 32-78. 

Silber, M. & Knobloch, E. 1989 Parametrically excited surface waves in square 
geometry Phys. Lett. A 137, 349-354. 

Simonelli, F. & Gollub, J.P. 1989 Surface wave mode interactions: effects of sym¬ 
metry and degeneracy J. Fluid Mech. 199, 471-494. 


